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ABSTRACT 



The analysis of multiple correlated two-dimensional (2-D) random signals or mul- 
tichannel 2-D signals is described. The emphasis is on estimation (linear prediction) 
and modeling of the 2-D random signals. Applications to spectrum analysis and image 
processing are considered. 
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I. Introduction 



This report deals with the analysis of multichannel two-dimensional (2-D) random 
signals. The signals are two-dimensional in that the signal is a function of two inde- 
pendent variables {i%i , n 2 ). We say the signals are multichannel since there are possibly 
several correlated 2-D signals, perhaps received on different communication channels 
that are processed and analyzed together. 

One example of a multichannel 2-D signal is the set of images received from a 
satellite multispectral scanner. The satellite forms several images of the ground in pixel 
registration but corresponding to different frequency bands in the visible and infrared. 
One can have four, seven, or even more such channels comprising the multichannel 
data. 

A more common example of a multichannel 2-D signal is a color image. The signal 
is inherently two-dimensional but consists of three registered components representing 
red, green and blue intensities or other tristimulus values such as those used in the 
NTSC video standard. 

Other examples of multichannel 2-D signals can be found in linear arrays when 
the sensor measurements are inherently separated along certain channels (for exam- 
ple in-phase and quadrature, principal and orthogonal returns of a radar, and in the 
measurement of dual polarized radar cross section as a function of space and time for 
extended targets. In these examples one can easily imagine multichannel signals that 
are higher than two- dimensional. 

Although the image processing community has long dealt with multiple correlated 
images, [see e.g. Ref. l], this community has for the most part not approached the 
analysis from a signal processing point of view. A notable exception is the work by 
Hunt and Kubler [2] who seem to have been the first to apply signal-processing methods 
to multichannel image restoration. The signal processing community on the other hand 
has dealt with the analysis of multidimensional signals but has also largely ignored the 
analysis of multiple multidimensional signals. While multichannel 2-D signals can be 
regarded as a special case of 3-D signals, multichannel signals have their own special 
properties that makes their analysis distinct and in many cases more tractible than the 
analysis of 3-D signals in general. Therefore we have chosen to concentrate on these 
special signals and study their properties explicitly. 

This report is primarily concerned with random multichannel 2-D signals and their 
applications. These signals will have properties sometimes similar to multichannel 
1-D signals and sometimes more like 2-D signals. We will be concerned here mostly 
with topics related to filtering and estimation of these signals; and in particular linear 
prediction and its applications. 

The remainder of this report is organized as follows. Chapter II deals with the rep- 
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reservation of multichannel. 2-D signals and their statistical characterization. Chapter 
III focuses on linear prediction and autoregressive (AR) modeling. This forms a basis 
for the remainder of the material. Chapter IV discusses spectrum estimation. Model- 
based methods for estimating the entire spectral matrix are presented. The estimate 
includes the 2-D autospectrum for each channel and the magnitude and phase of the 
cross- spectra. Chapter V describes applications of the theory to image processing. 
Chapter IV provides a summary and conclusion. Areas for future work are cited there. 
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II. Statistical Characterization of Multichannel 2-D Random Processes 
2.1 Representation of Multichannel 2-D Signals and Systems 
2.1.1 Signal Domain Representation 

Discrete multichannel 2-D signals (or sequences) will be represented by a vec- 
tor quantity x(n!,n 2 ) as illustrated in Fig. 2.1. Each component x K (n 1 ,n 2 ) of the 
vector represents a 2-D signal existing in one of the channels. Note that the vector- 
valued nature of the function relates to the fact that the signal is multichannel while 
the vector nature of the argument n = (nj,n 2 ) relates to the fact that the signal is 
multidimensional * 

The processing of these signals is analogous to that for other discrete signals. We 
will be primarily concerned with linear shift-invariant (LSI) operations that can be 
represented by 2-D vector difference equations. In particular an M-Channel 2-D signal 
x(ni , n 2 ) is transformed to another M- channel 2-D signal by an LSI system represented 
by the difference equation 

y(ni,n 2 ) = — A Xli .y(n 1 -i 1 ,n 2 -i 2 )+ B tlll x(n 1 - l x , n 2 - l 2 ) 

(£i .£ 3 ) e <3 

(« 1 .» 3 )*( 0 , 0 ) 

( 2 . 1 ) 

where the A ili2 and B lltl are M x M matrix coefficients and a and (5 define their 
(fixed) regions of support. The system will be recursively computable if there exists an 
ordering for processing of the points such that values of the signal y needed to compute 
the signal’s current value at the point (n 1? n 2 ) are always available. Whether such an 
ordering exists depends on the shape of the output region a. 

The output for a multichannel 2-D system can be equivalently represented by the 
2-D convolution operation in either of the forms 

CO , 00 

y(ni,n 2 )— H (ti , , —li , n 2 — £ 2 ) (2.2a) 

£ j ,£ 3 = — * 00,-00 

CO , oo 

y(ni,n 2 )= ^ H (n x - , n 2 - £ 2 )x(£i , t 2 ) (2.26) 

£ 1 £ 3 — — CO , — OO 

where H (•, •) is a matrix function representing the 2-D multichannel impulse response. 
(Each term if(£ l5 £ 2 ) in Eq. (2.2) is an M x M matrix.) A 2-D multichannel system 



* Although it is tempting to use vector notation for the arguments, this tends to 
hide the essential 2-D nature of the signals. Therefore we use the longer but more 
explicit notation involving the two indices n x and n 2 . 
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Channel m 




x 1 (n ) ,n 2 ) 
X^, n g ) 

X 3 (n l ,n 2 ) 



Figure 2.1 Discrete Multichannel 2-D Signal. 
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will be called a Finite Impulse Response (FIR) system if the support of is finite 

and an Infinite Impulse Response (HR) system otherwise. Clearly for an FIR system 
all of the A ili7 in Eq. (2.1) are zero and H{t l ,t 2 ) = B lll2 for (£i,£ 2 ) € /?. 

For purposes of this report, we will be dealing with statistical properties of random 
signals. Thus the mean of a signal x(n!,n 2 ) is a vector quantity defined by 

m* («i , n 2 ) = £[x(n! , n 2 )] (2.3) 



and the correlation and covariance are M x M matrix functions defined (respectively) 
by 

R(n l ,n 2 ; m 1 ,m 2 ) = E [x(n 1 ,n 2 )x T (m 1 ,m 2 )] (2.4a) 



C(n 1 ,n 2 ;m 1 ,m 2 ) = E (x(n!,n 2 ) - mjtinnj)) (x(m 1 ,m 2 ) - m I (m 1 ,m 2 )) J 



(2.46) 



When the random process representing the signal is homogenous or stationary the mean 
is constant and the correlation and covariance are functions only of the vector distance 
between the points (n l5 n 2 ) and (m 1 ,m 2 ). In this case Eqs. (2.3) and (2.4) become 



and 



m* = E [x(n x , n 2 )] 



(2.5) 



R(ki , k 2 ) — E 
C(k l2 k 2 ) = E 



[x(n 1 ,n 2 )x r (n 1 -k ly n 2 - k 2 )] 
(x(nj , n 2 ) - m x ) (x(n x - k 1 , n 2 



- k 2 ) - m,) T 



(2.6a) 

(2.66) 



It can be seen from their definition that the correlation and covariance functions have 
the following symmetry properties 



R{ki,k 2 ) = R T (-k 1 ,-k 2 ) (2.7 a) 

C{k 1 ,k 2 ) = C T {-k 1 ,-k 2 ) (2.76) 



When a stationary random signal is transformed by a linear shift- invariant system, 
the statistical characteristics of the output can be expressed in terms of those of the 
- input by using Eqs. (2.1) and (2.2) and taking expectations. Thus, from (2.1) the 
mean satisfies 

( * x i *’ 3 ) € a ( l i » 1 3 ) £ P 

(*1 ( 0 , 0 ) 



or 



- i 



7~ , ^.1‘j 



E 



,('142)60 



2 



m T 






(2.8) 
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where 



I 



(2.9) 



The correlation function satisfies the pair of equations 

R y (k i,k 2 ) = — Y, A il{l R y (k 1 —i 1 ,k 2 —i 2 )+ Y, B ll(2 R xy (k 1 — , k 2 —i 2 ) 

(*li*2)€ Of (^11^2)6^ 

(*i.* 2 )^( 0 , 0) 

(2.10a) 

R yx {k 1 5 ^ 2 ) “ ” ^ ^ t 3 (^1 ? ^2 “~*2 ) “t* ^ ^ Rt\ iti Rx (^1 “^1 5 ^2 ”^2 ) 

(*11*2) 6 a (^11^3)6/? 

(*1 »*2 (0,0) 

(2.106) 

where the cross correlation is defined by 

R xy {ki,k 2 ) = R yx {-k 1 ,-^ 2 ) = E [(x(n 1? n 2 ) - m,)(y(n 1 ,n 2 ) -m y ) T ] (2.11) 
Alternate relations, in terms of the impulse response are, from (2.2) 



m y 




H{1 1 , 4 ) 



( 2 . 12 ) 



and 

00 , 00 00 , 00 

Ry(k 1 ,k 2 )= Y H(i 1 ,t 2 )R x (k 1 +p 1 -£ 1 ,k 2 +p 2 -l 2 )H T (pi,p 3 ) 

£ 1 , £ 2 — — 00 , — 00 pi,p 2 = -oo,-oo 

(2.13) 

The relation between the correlation and the cross correlation can also be expressed in 
terms of the impulse response as 



OO , OO 

Ryx{ki,k 2 )= Y H(i 1 ,t 2 )R x (k 1 - £i,fca - la) (2.14 a) 

(£i,£ 3 )= — 00,— 00 

and 

OO ,00 

Ry {k 15 ^ 2 )“ ^(^l)4)^iy (^1 ” ^ 15^2 “ 4) (2.146) 

£ l : £ 2 — — OO , — OO 

which forms a connection between Eqs. (2.10) and (2.13). Relations analogous to Eqs. 
(2.10), (2.11), (2.13), and (2.14) clearly apply to the covariance and are easily derived. 

2.1.2 Frequency Domain Representation 

When a (deterministic) signal of the form 

x(n l5 n 2 ) = cz^z^ (2.15) 
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(2.16) 



is applied to an LSI system, the output (from Eq. 2.2) is given by 



y(n l5 n 2 ) = H z (z 1 , z 2 )cz^ z ” 5 



where 

co , co 

H z (z x , z 2 ) = ^ H{h,l 2 )z- tl z; 1 ’ (2.17) 

, £2 = “ OO , — OO 

is the z transform of the impulse response and is called the system function . If one 
applies Eqs. (2.15) and (2.16) to Eq. (2.1) it can be seen that 



H z ( Zi , z 2 ) 



E 



Aj { z, 
1 



. ( * 1 , * 3 ) € a 



y 1 

j 



E 

Ji.hep 






-h 




(2.18) 



The frequency response of the system is the system function evaluated on the two unit 
circles 

H, u {u 1 ,oj 2 )=H z (2.19) 

Stability of the 2-D multichannel system corresponds to convergence of H z on the unit 
circles. Stability thus implies that the frequency response exists for all values of uq and 
U 2 • 



Random signals are characterized in the frequency domain by the power spectral 
density matrix. This is defined as the 2-D Fourier transform of the correlation function 



OO . CO 

■Muq,^) = Y, R x {k l ,k 2 )e- j “' k 'e~ j “'- k ' (2.20) 

/c ! ,/c 3 = — 00,— CO 



The spectral matrix is positive semi-definite and Hermitian. Its off-diagonal terms 
represent cross spectra between components of the signal and need not be real. For the 
case of a two-channel 2-D signal the matrix can be written as 

or,, , , \ _ r-S'n(wi,w 2 ) S 12 (uq,w 2 ) 

^ x 5 ^2 J c t ( \ O ( \ 

_^21 ( W 1 5 W 2 j *522(^1 5^2) 

The terms S X1 and S 22 are real and nonnegative and represent the 2-D power spectrum 
or autospectrum for each channel while the terms S 12 and S 21 represent cross spectra. 
Because of the Hermitian symmetry the cross spectra have identical magnitude and 
opposite phase. The normalized quantity 



(2.21) 



/C 2 (uq ,U> 2 ) = 



l‘S'12 ( w i > w 2 ) | 2 

■Si ! (wj , u; 2 ) S 22 (uq , u > 2 ) 



(2.22) 



is the 2-D magnitude squared coherency (MSC) and is frequently used instead of |S' 12 | 
to represent the magnitude of the 2-D cross spectrum. Further treatment of the MSC 
is given in Appendix A. 
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A cross - spectral matrix S xy (u 1 ,u 2 ) can also be defined as the Fourier transform 
of the cross correlation function 

oo ,oo 

S„(w.,w») = £ (2.23) 

/c J. , /c 3 = — CO , — OO 

Its components represent cross spectra between components of the signal x and com- 
ponents of y. This matrix is neither positive definite nor Hermitian. However it shares 
a symmetry property with S yx in that 

= (2.24) 

When a signal is transformed by an LSI system, frequency domain expressions for 
the output power spectrum and the cross spectrum between input and output can be 
obtained. These are derived by taking Fourier transforms of Eqs. (2.13) and (2.14). 
The output power spectral matrix has the form 

S.J (u>. , Un ) = H. v (u. , V; ) S x (uq , Un ) H* T (l- 1 , Wo ) 

while the cross spectral matrices relate to S x and S y through the expression 

Syx ( w i ) ) — H w (cuj ,u> 2 )S X (uq , ui 2 ) ( 2 . 26 a) 

s y (wj , u 2 ) = H w (wj , u 2 )S xy (a?! , u> 2 ) (2.266) 

and the symmetry relation (2.24). 

2.2 Vector Representation of Multichannel 2-D Signals 
2.2.1 Signal Vector Representation 



(2.25) 



Frequently, a multichannel signal is defined over some finite rectangular region of 
support. That is, the signal is considered only for 0 < < N l , and 0 < n 2 < N 2 . 

In this case it can be useful to represent the signal as an N X N 2 M - dimensional vector 
of the signal values. Two such representations are most useful. In the first, the vector 
valued signals at points ( rii,n 2 ) are lexicographically organized first by row and then 
by column into the larger vector. That is, the multichannel 2.-D signal is represented 

by 

(2.27a) 




where 

‘ x(n l5 0) ‘ 

x(n a , 1) 

— " i = ; 

.x(n! , N 2 - 1). 



(2.276) 
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This representation will be called index ordering . 



A second type of representation organizes signal component within each channel 
by rows and columns and then stacks the components. This representation has the 
form 



where 



x = 



X = 



-i 

^0 



— Ni - 1 “ J 



X 1 



,M 



m — 1,2, ... ,m 



where 



x = 



K,0) 

(^1,1) 



n i — 0 | 1 > • • • > Ni 1 



(2.28a) 



(2.286) 



(2.28c) 



■Im (^1,-^2 - 1) 

This representation will be called component ordering . 

Index and component orderings are related by permutation transformations 



x' = Px 
x = P r x' 



(2.29a) 

(2.296) 



The form of the permutation matrix is illustrated in Fig. 2.2 for N 1 =3 , N 2 = 4, and 

M = 2. 



2.2.2 Statistical Characterization of Signal Vectors 

The signal vectors previously defined can be characterized by mean vectors and 
correlation or covariance matrices . These quantities are defined for index-ordered vec- 
tors by 



m = P[x] (2.30a) 

R = F[xx r ] (2.306) 

K = P[(x — m) (x — m) r ] (2.30c) 

and for component-ordered vectors by similar expressions with primed variables. From 
Eqs (2.29) and (2.30) it can be seen that the index-ordered and component-ordered 
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X 



Figure 2.2 Permutation matrix for Nj = 3, N 2 = 4. M = 2 . 



statistics are related by 



(2.31a) 

(2.316) 

(2.31c) 



m' = Pm 
R' = PRP r 
K' = PKP t 



Note that for a stationary random process the vector m will consist of N x iV 2 M- 
dimensional subvectors equal to the signal mean (see Eq. 2.5). The vector m' will 
consist of MiY x iV 2 -dimensional subvectors. Each subvector has elements that are all 
identical and equal to the mean of the signal in the corresponding channel. 

The correlation and covariance matrices for a stationary 2-D random process have 
specific structure and symmetry at their various levels of partitioning. These properties 
are illustrated in Figs. 2.3 and 2.4. An example showing detailed structure of these 
matrices is given in Appendix B. 

2.2.3 Linear Transformation of Signal Vectors 

Linear transformations on signal vectors can be represented by matrix equation 

y = Ax (2.32) 

We will consider only index-ordered transformations here since index ordering is more 
extensively used in this report. The equivalent matrix for component-ordered transfor- 
mation is given by 

A' = PAP t (2.33) 

General linear transformations such as this may operate upon the rows, columns, or 
channel dimensions of the signal. 



Linear transformations corresponding to a LSI filtering operations are represented 
by a matrix of the form 



H = 



H(0) H(— l) . . . Ht-iVi+l) 

H(l) H(0) . . . H(— iW+2) 



_H(iV a -l) H(iW-2) 



where 



HK) = 



iL(ni ,0) 

H [n l , l) 



1 ) 

H{n i,0) 



lH(ni,N 2 — l) H{n 1 ,N 2 -2) 



H(0) 



H{n lt -N,+ 1) 

H (n : , — iV 2 + 2) 

H(n ! , 0 ) 



(2.34a) 



(2.346) 
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Signal Vector 



Xo 




x(n,.0) 


*1 




x(n,,l) 




where x n = 






— r *i 




X Yl -l 




x(n,.;Y,-l) 



Correlation Matrix 



RjO) R(l) 

R(-l) R(0) 



R=£ = 



R(A'i-l) 
. . R(.V 1 -2) 



[BLOCK TOEPLLTZ WITH 
.\\MxX 2 M BLOCKS) 



R( — jV, — 1) R(-A r ,-2) R(0) 



where 



R(k.O) R(k.l) . R[k..\\-1) 

R(k-.~ 1) R[k,0) . R(k..\\- 2) 



R(*) = £ x n x n r _* 



[BLOCK TOEPLITZ WITH 
M x M BLOCKS) 



R[k.-.\\- 1) R[k- A'j-2) R(k.O) 



where 



R[k l .k i ) = R T [-k l -k,) = Ex[n l ,n,)x T (r h -k i .n,-k 2 )\ [HOT TOEPLITZ ) 

Fig. 2.3 Correlation matrix for index ordering. 
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Signal Vector 



X 1 




-y ™ 
±0 




i*(ni-0) 


o 

X" 




XT 








where x m = 




where x™ = 

- n i 




X M 




v rn 

X.N'-l 










— l 







Correlation Matrix 



R.i 


R 12 


R 


R 2 J 


Ro o 


R 



R ' = £ x 'x T 



(.V, A'.xAVVj BLOCKS ) 





• ■ R-m.w 



where 



0/ k r>l k Dt k 

D Ik nl k r>l k 

n_i n 0 ... rc Vi _o 



R; A = £ x'x* 7 = 



(BLOCK TOEPLITZ 
WITH S 2 > N, BLOCKS) 



pi k 



R 



/ A 
— iV. 



R 



i k 

o 



where 



R 



t k 

n, 



-f -y-' V* r 
- Cj x„x r> _„ i 



[TOEPLITZ) 



Fig. 2.4 Correlation matrix for component ordering. 
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(The ~ notation has special meaning that will be explained shortly.) When the trans- 
formation corresponds to an FIR filter many of the blocks and subblocks in Eq. (2.34) 
will be zero; the specific arrangement will depend on the shape of the support region 
for the impulse response. 

For a system described by a linear difference equation (2.1), one can write 



Ay = Bx (2.35) 



where the matrices A and B are defined in terms of the coefficient matrices and 

Bj analogously to Eqs. (2.34). The output y can then be expressed directly as 



y = 




(2.36) 



This form of the analysis implicitly assumes zero initial conditions for the difference 
equation. 

Occasionally it is necessary to convert between one form of channel representation 
for the 2-D multichannel signals and another without performing any specific filtering 
or signal processing. This would occur for example if the signal x(n.! , n. 2 ) represented 
a color image with red, green, and blue components and we wanted to convert the 
signal to NTSC form. If the new signals are defined in terms of the old signals by the 
transformation 

x(n lt n 2 ) = Tx(ni,n 2 ) (2.37) 

then the index-ordered signal vector would be transformed according to 



x = Tx (2.38) 

where T is a N l N 2 M by N X N 2 M block diagonal matrix. This matrix would contain 
N ! N 2 nonzero blocks all equal to the signal transformation T. 

2.2.4 Reversal Notation 

In the analysis of 2-D signals using vector representation, it is sometimes necessary 
to deal with vectors and matrices whose components are ordered corresponding to 
decreasing values of the signal index rather than increasing values as in Eq. (2.27). 
While a reordering of the signal values can be represented by a permutation matrix 
operation, it is much clearer to develop explicit notation. The notation to be used 
will apply only to index-ordered vectors and matrices. Component-ordered quantities 
are less frequently used and do not lend themselves readily to the reorderings that are 
considered here. 
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The symbol ~ above a vector will correspond to a vector whose first-level partitions 
have been organized according to decreasing value of the index. 



x — 



X Wl -i 

X /^-2 



*0 



(2.39) 



where the x„ are defined by Eq. (2.27b). We call this a first index reversal of the 
vector. The symbol ~ below a vector will be used to denote a vector whose second 
level partitions have been organized according to decreasing index values. 



*1 



LXjV.-l - 



(2.40) 



where 



x = 



x(n! , N, — 1) 
x(ni ,N 2 — 2) 



(2.41) 



x(jt, ,0) 

We call this a second index reversal of the vector and is occasionally useful. Most 
frequently it will be required to use the doubly reversed vector which is denoted by a 
double~over the vector: 



x = 



2^-1 

2^-3 



*0 



(2.42) 



and denotes a vector where both the main blocks and the second level blocks have been 
reversed. 



Since transformation of a reversed vector usually implies generation of a new vector 
with corresponding order, it is most reasonable to define reversals of a matrix as a 
reordering about both its rows and columns. Thus a matrix A represents a matrix A 
whose first level blpcks or partitions have been reversed along both rows and columns. 

Matrices and A are defined similarly with reversals about rows and columns of the 
inner blocks. Note also that reversal about both rows and columns of a matrix is 
equivalent to transposition about both the main and reverse diagonals. This preserves 
certain symmetry such as block Toeplitz structure when it exists. 



If reversal notation is used, the LSI transformation of Eq. (2.2) can be represented 



as 



y = Hx 



(2.43) 
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where the matrix H is now defined by 



with 



H(nO = 





H(0) 


H(l) 


... H(iV 1 - 


H = 


H(— 1) 


H(0) 


... H(A r 1 — 




.H(— iVi + 1) H(— IVi + 2) 


H(0) 


■ H(ny, 0) 


H{n lt 1) 


... H{n lt N 2 


1) 




... H (riy , N 2 


-H(rii,N 2 — 1) 


H{ny,N 2 -2) 


H (ni , 



(2.44a) 



(2.446) 



Reversal notation will also prove to be useful in considering signal models in Chapters 
III and IV. 

2.3 Separable Multichannel 2-D Signals and Transformations 
2.3.1 Separable Signals 

The direct product of an L x Q matrix B and an N x P matrix A is an NL x PQ 
matrix 

B cl 1 1 B cl 1 2 ... Ba y p 

Ba 2l Ba 22 ... Ba 2 p 



B ® A — 



(2-45) 



, Ba^ i Ba^j 2 ... Ba ^ p _ 

In general B <g) A is not equal to A 0 B. A number of other properties are given in 
Table 2.1 and in Appendix C . 



Table 2.1 - Properties of Direct Product 



(B ® A) (D ® C) 


=BD ® AC 


{B®A) T 


=B T 0 A T 


(BOA)' 1 


=B ~ 1 0 A -1 



Consider for the moment a single channel 2-D random signal that has the separable 

form 

x(ni,n 2 ) = x il (n 1 ) • x B (n 2 ) (2.46) 

Then if x A and x B are the vector representations of the signals (n^) and x B (n 2 ), a 
vector representation of the 2-D signal is given by 

X = x B ®X A ( 2 - 47 ) 
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Now if (nj ) and x B ( n 2 ) are independent random processes with correlation matrices 
K a and R B the correlation matrix for the 2-D signal vector is 



E [xx T ] = E (x B 0 x A ) (x B 0 x A )' 

= E [(x b Xb)] ® E [(x A x^)] 



(2.48) 



— R B 0 R>i 

A similar direct product relation holds for the mean vector and the covariance matrix. 



Multichannel 2-D random processes can be separable along the index (nx,n 2 ) 
directions, between channels, or both. Various cases of separable signals and their cor- 
relation matrices are given in Table 2.2 assuming index ordering. The specific structure 
of the matrices can be envisioned by comparisons with the matrices in Fig. 2.3. 



Table 2.2 - Separable Forms for Multichannel 2-D Signals 



Form of Signal 


Description 


Form of Correlation 


x(n.j , n 2 ) = x A (n x ) -x B (n 2 ) 


Product of a 1-D 
single-channel 
signal and a 1-D 
multichannel signal 




x(nj,n 2 ) = x(n 1 ,n 2 ) • c 


Product of a 2-D 
single channel signal 
and a M-dimensional 
random vector 


R^ = R c 0 R x 


x(n.j , n 2 ) = x A (nj • x D (n 2 ) -c 


Product of two 1- D 
single channel signals 
and a random vector 


R^ = R c 0 R b 0 R^ 



Separable structure for the signals when it exists leads to simplifications in the 
analysis. For example 2-D problems can be reduced to a pair of 1-D problems. The 
implications of separability for linear prediction are described in Chapter III. 

2.3.2 Separable transformations 

If a linear transformation for a multichannel 2-D signal is completely separable it 
can be represented as the direct product of three matrices 

T = W 0 V 0 U (2.49) 
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where U is a Ni x Ni matrix representing transformation along the n x direction, V is an 
N 2 x N 2 matrix representing transformation along the n 2 direction and W is an M x M 
matrix representing transformation between channels. If the signal is correspondingly 
separable as 

x = c®x B ®x^ (2.50) 

then the processing of the signal is greatly simplified since 

Tx = ( Wc ) ® (Fx B ) ® (Ux A ) (2.51) 

However separable transformations are advantageous even if the signal is not separable. 
Common examples of separable transformations are the DFT, Hadamard and Walsh 
transforms, and the singular value decomposition. 

A convenient relation exists between index-ordered and component-ordered sepa- 
rable linear transformations. In particular, if the index-ordered transformation is given 
by (2.49) then the component-ordered transformation is given by 

T = V ®U ®W (2.52) 



A somewhat more concise representation for separable transformations exists. This 
involves first representing the multichannel 2-D signal in yet another form, as a N X M x 
N 2 matrix 



r*i i 




lx M ] 



(2.53) 



where X m is a x N 2 submatrix whose elements are the signal x m [n l , n 2 ) in channel 
m. It can be verified that with this representation, a separable transformation of the 
form 

y = (W ® F ® 17) x (2.54) 

leads to a matrix representation Y with blocks 



M 

Y m =J2^meUX e V T 
£= 1 



m = 1, 2, . . . , M 



(2.55) 



where w ml are the elements of W . This can be written in a single matrix equation as 



Y = (/ ® W) ■ UXV T 



(2.56) 
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III. Linear Prediction for Multichannel 2-D Random Processes 

3.1 Linear Prediction 



3.1.1 Problem Statement 

Let x(n 1 ,n 2 ) represent a zero-mean stationary multichannel 2-D random signal. The 
linear prediction problem is concerned with forming a linear estimate of the signal at the 
point (n l5 n 2 ) from other values of the signal in a region a. Specifically we form the 
estimate 

x(n l5 n 2 ) = - ^ A? i ^{n 1 — h,n 2 — 1 2 ) (3.1) 

)€<* 

(«i,»j)*(0,0) 

and the prediction error is defined as 

e(n 1 ,n 2 ) = x(n l5 n 2 ) - x(n l5 n 2 ) (3.2) 

The matrix coefficients Aj ; are chosen to minimize the mean- squared error E ' e ( n. j. . ) | 3 . 

The error is generated from the data by the FIR prediction error filter. From Eqs. 
(3.1) and (3.2) the difference equation for the filter is 

e(n l5 n 2 )= A J 1 < i x(n 1 - t L ,n 2 - t 2 ) (3.3) 

( *1 i * 7 ) 6 Of 

with .Too = I ■ Since the filter is FIR, the terms Aj^ . also represent the impulse response 
and a is the impulse response region of support (see Eq. 2.2a). The correlation function 
of the error evaluated at lag (0,0) is the M x M prediction error covariance matrix E,. 
That is 

= E [£(n 1 ,n 2 )e r (n l5 n 2 )] = R e ( 0,0) (3.4) 

The mean-squared error £ is given by 

£ = E [|c(n 1 ,n 2 )| 2 ] = E [e T (n x , n 2 )e(n 1 , n 2 )] = trE e (3.5) 

These quantities are important for the analysis that follows. 

3.1.2 Normal Equations 

The prediction error filter coefficients and the prediction error covariance are formed 
by solving a set of Normal equations . These linear equations follow directly from the 
orthogonality principle [4] which states that the error must be orthogonal to the signal 
values used in prediction. 
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(3.6a) 



E [x(n x — i'x , n 2 - t 2 )e r (^ 1 , n 2 )] = [O] 

(j'i,z 2 )€a:, (*i,t 2 ) ^ (0,0) 

It further follows from Eqs. (3.3), (3.4) and (3.6a) 

E [x(n 1 ,n 2 )e r (n 1 ,n 2 )] = E [e(n x , n 2 )e r (n x , n 2 )] = E £ (3.66) 

These last two equations can thus be written more concisely as 

E [x(n x — fi,n 2 — f 2 )e T (n l5 n 2 )] = E e 6(i 1 )<5(i 2 ) 

(i'i , f 2 ) G Q! (3. / ) 



3. 1.2.1 Rectangular Support 



The Normal equations will be developed by using an index-ordered vector representa- 
tion for the signal. At first let a be the rectangular region shown in Fig. 3.1. The size of a is 
P l xP 2 points and Lx and L 2 can be any values in the range —Pi < L x <0, — P 2 < L 2 < 0. 
Then an index-ordered representation for the signal data that occurs in Eq. (3.3) is needed. 

To make the analysis clear, consider the case where L y = L 2 = 0. This is the quarter 
plane or first quadrant predictor. The signal data needed in Eq. (3.3) is in the range 
rix — Pi + 1 to n 1 and n 2 — P 2 + 1 to n 2 (see Fig. 3.2). Form a vector x ni „ 2 as follows 



w n i n 2 



-ni-P i + l 

— n i — P i + 2 



x — P x + i 



x[n l - Pi + t, n 2 — P 2 + 1) 
x(ni - Pi + i, n 2 - P 2 + 2) 



x(ni - Pi + t , n 2 ) 
Then Eq. (3.3) can be written in matrix form as 

where 



n 2 ) 


= A r x „ 

— n i n 




r a<°> i 


— 


A (1) 








-4t0 


— 


-4a 




-■4>,f 3 -i - 



(3.8a) 



(3.86) 



(3.9) 



(3.10a) 



(3.106) 
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Figure 3.1 A general rectangular region of support ac 
for the prediction error filter. 
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V’ 

p ,-' 



I a" X(n-i 



i ,= 0 V2 

l 2~ 0 



Figure 3.2 First quadront predictor. 
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The Normal equations then follow directly. From Eqs. (3.7) through (3.9) we have 











E 


x rllrt5 € T (n 1 ,n 2 ) 


= E 





or 



where 



= E 



x„ „ x„ „ 

— n irii — n j rio 



A = S 



(3.11) 



RA = S 



(3.12) 





(3.13a) 



(3.136) 



and the partitioning in Eqs. (3.13) corresponds to that of Eqs. (3.8) and (3.10). Eq. (3.12) 
represents the Normal equations. The results for a first quadrant multichannel predictor 
are summarized in Fig. 3.3. 



Observe that the correlation matrix that occurs when the Normal equations are writ- 
ten in the standard form (3.12) is the reverse matrix R and not R (compare Figs. 2.3 
and 3.3). This point is not usually appreciated or discussed in the literature on linear 

prediction. For the single channel case R and R are both the same; for our multichannel 
case it is important to make the distinction. 



When L x and L 2 of Fig. 3.1 are not equal to zero a similar analysis again leads to 

Eq. (3.12). The correlation matrix R is identical to the one for the previous case but the 
matrix of filter parameters is defined by 



A = 



A (Ll) 

L i + P i - 1 ) 

A\L 3 



A (i) = 



^t,0 






+ P 2 — 1 _ 



(3.14a) 



(3.146) 
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Prediction Equation 



r 1 — 1 / 3 — i 

x(n L n 2 ) = 2 2 (i„t,H0.0) 

Ii = 0 t 2 = 0 



Normal Equations 



R(0) R(-l) 

R(l) R(0) 


• R(-P,-1) 

R(-P,-2) 




O — 

.<< ■ 




0 VL 

0 


R(P.-1) R(P|- 2 ) • 


R(0) 




A i p .-" 




0 



where 





R{k,0) 


R(jt.-l) 


. R(k,-P 2 -i) 








S e 




fl(M) 


R(k.O) 


. R(k,-P 2 ^ 2) 




-Vl 




0 


R(k)=R r (-k)= 








A<*> = 




s (0| = 








R(k.P,-2) . 


. . R{k. 0) 




Ak.p t -i 




0 



where 

( Atj . ^3) = T (-k l ,-k 2 )- E\x(n i 1 n 2 )'K T ( A 0 0 =I E f =E (x-x)(x-x) r ] 



Fig. 3.3 Equations of linear prediction for first quadrant support. 
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The matrix S is similar to Eq. (3.13) but the non-zero block appears in a position 
corresponding to the location of the coefficient A 0 o in the matrix A. 

3. 1.2. 2 Special Forms 

The Normal equations for rectangular support can be viewed in the following way. 

We have a set of Pi x P 2 points as shown in Fig. 3.4(a). A correlation matrix R is formed 
for the P X P 2 points. The terms of the correlation function appearing in this matrix are 
shown in Fig. 3.4(b). The Normal equations for predicting various points in the array all 

involve this same matrix R. They differ only in the indexing of the filter coefficients and 
in the position of the terms A 0 o and in the arrays A and S. 

The filters for predicting the upper right and upper left points in a rectangular region 
of support are called the first quadrant and second quadrant filters. The first quadrant 
filter was discussed in detail in the previous subsection. By using properties of the reversal 
operation the second quadrant Normal equations can be put in the form of Fig. 3.3. When 
this is done it can be observed that the equations although similar have their second level 
blocks transposed. Thus the filter coefficients for the first and second quadrant predictors 
are not the same. The third and forth quadrant filters for predicting the bottom left and 
bottom right points are also distinct. Their Normal equations differ from the Normal 
equations of the first and second quadrant filters in that the innermost blocks R(ki,k 2 ) 
are replaced by their transposes. This is a difference from the single channel 2-D case 
where the first and third and the second and fourth quadrant filters are identical. 

A very important form of predictor is the nonsymmetric half plane filter. Although this 
filter clearly has non-rectangular support the form of the Normal equations can be derived 
from those already considered. The linear prediction problem is depicted in Fig. 3.5. One 
can think of the NSHP support region as a rectangular region where |L 2 | points (marked 
with x’s in Fig. 3.5(a) are missing. The circled point is the one predicted. The Normal 
equations for the NSHP can be obtained by starting with the equations for rectangular 
support. One then simply drops out filter coefficients corresponding to the “missing” points 
and eliminates corresponding rows and columns of the correlation matrix. The required 
terms of the correlation function are shown in Fig. 3.5(c). The detailed form of the Normal 
equations is given in Fig. 3.6. 

A companion NSHP prediction problem is shown in Fig. 3.5(b). While in Fig. 3.5(a) 
data is processed from lower left to upper right, in Fig. 3.5(b) data is processed in the 
opposite direction. We refer to Fig. 3.5(a) as forward prediction and Fig. 3.5(b) as 
backward prediction. The Normal equations for the backward problem can be obtained in 
the same manner discussed for the forward problem. If the equations are permuted and 
put in the form of Fig. 3.6, they will differ only in that the inner blocks R(ki,k 2 ) appear 
transposed. 
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(a) 



P -1 



R(k ] , k 2 ) 



-P. + 



V 1 



- p i + l 



(b) 



Figure 3.4 Linear prediction with rectangular support. 
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Note : l_ 2 is a negative number 

Figure 3.5 Linear prediction with NSHP support. 



Prediction Equation 



u-P >- 1 _ R.-i 

x(n 1 n 2 )= £ -•l 0 \ J x(n 1 .n 2 -i 2 )- E E ^ x ( ” l ~ 2 1 * ~ * 2) 

tj = 0 ! 1 = 1 U = ^: 



Normal Equations 



R ( 0 ) R'(-i) 

R'U) R(°) 


. R'(-P, + l) 

. . R(-P,+2) 




A d°) 
A (» 




S ' (0) 

0 


R'(P,- 1) RIP,- 2) . 


. R(0) 




A-" 1 '; 




0 



where 





R{h-L,) 
R(k,-L 2 - 1) 


R(k.-Ln-l) 

R[k,-L 2 ) 


R(k,~ L 2 - P 2 -rl) 
. R(k-L 2 -P 2 ~ 2) 




A k.L 2 

A k,L 2 +l 


R'(Jfc) = R' r (-*) = 


R{k.-L 2 -P 2 - 1) 


R(k.-L 2 ^P 2 - 2) . 


R(E0) 


A (*) = 


A k,L 2 ~P 2 -l 



and 



R ( 0 ) = 


R(0,0) 

R(0,1) 


R(0.-1) 

R(0.0) 


. . R(0.-L 2 -P 2 -l) 

. . R(0,-L 2 -L 2 -P 2 ^2) 


A'(°' = 


0 0 

0 


S'(°> = 


0 




R{Q.L 2 -P 2 -l) 


R(0.L 2 - P 2 -2) . 


R(0.0) 








0 



and where all other quantities are defined as in Fig. 3.3. 



Fig. 3.6 Equations of linear prediction for NSHP support (forward predictor). 
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3.2 Autoregressive Models 



If the error process e(n.! , n 2 ) produced by linear prediction were white then one could 
obtain the original random process by inverse filtering. Specifically, by replacing e(n. l5 n 2 ) 
by a white noise w and solving Eq. (3.3) for x(n.!,n. 2 ) we have 

x(n 1 ,n 2 ) = - ^ Af i<3 x(^i - i\,ri 2 - h) + w(n x ,n 2 ) (3.15) 

(»' l .<3 )€ a 
(»i, *3)^(0, 0) 

This is the multichannel 2-D autoregressive (AR) model. The white noise has a covariance 

E w = E [w(n 1 ,n 2 )w T (n 1 ,n 2 )] (3.16) 

Fig. 3.7 illustrates the relation between linear prediction and autoregressive modeling. 
While linear prediction uses an FIR filter, AR modeling employs' a recursive or HR filter. 

Since the filter coefficients in an AR model satisfy Normal equations identical to those 
for linear prediction, there is a temptation to equate the two concepts. The distinction 
lies in the fact that linear prediction does not usually produce white noise. When one uses 
an AR model to represent an arbitrary random process and sets up Normal equations to 
solve for the filter coefficients, one is essentially ignoring the distinction. The justification 
for the white noise AR model can lie only in the belief that an AR model of some order is 
a close approximation to the true process. For a model with NSHP support it is true that 
the AR model becomes a close approximation to any given random process if the order 
becomes sufficiently large. This argument follows by generalization of results in Refs. 5- 
and 6 to the multichannel case. 



3.3 Linear Prediction for Separable Problems 

When the 2-D multichannel random process is separable, various simplifications result 
in the linear prediction problem. Because the correlation matrix is separable (see Section 
2.3), the Normal equations (3.12) for linear prediction problems with rectangular support 
reduce to a set of lower-dimensional Normal equations. The results for a predictor with 
support in the first quadrant are summarized here. 



When the correlation within channels is separable from the correlation between chan- 
nels the covariance matrix R in (3.12) can be represented as the direct product of a 

between-channel covariance matrix R c and a within-channel covariance matrix R x . In 
this case we can postulate a separable form for the matrices A and S and show that they 
provide a solution to the original Normal equations. In particular, we can write 



(Rc ® R* ) (I m ® a* ) 




(3.17) 



where I M is the M x M identity matrix, a* is an N 2 -dimensional vector whose first 
component is 1, and is the Ni N 2 -dimensional vector (l, 0, 0 . . . , 0) T . Then (3.17) 

implies the two relations 



R*a x = a\t Nx 



N, 



(3.18a) 
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(a) Linear prediction 
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(b) AR modeling 



Figure 3.7 Linear prediction and A R modeling. 
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E ( =o 2 z R c (3.186) 

Since the solution of (3.12) is unique, this shows that it can be found from (3.18). The 
form of (3.17) indicates that each channel has the same filter coefficients a x and that there 
are no between-channel terms. The parameters a x and a\ are obtained by solving the 
Normal equations (3.18a) pertaining to a single channel 2-D problem. Then from (3.18b), 
the prediction error covariance for the 2-D multichannel problem is a scaled version of 
the between-channel covariance. The scale factor is the single channel prediction error 
variance. 

Similar results can be obtained if the correlation matrix is separable along the n l and 
n 2 directions. In this case the correlation matrix can be represented by the direct product 
of a matrix R B which relates to a 1-D multichannel process along the n 2 direction and 
a matrix representing a 1-D single channel process along the n l direction (see Table 
2.2). The postulated solution is of the form 



/ i \ 

(R b ®R*)(A b ®a*) = ( — S ,0) 1 

\ a \ J 



& 



-.v. 



(3.19) 



which implies 



R B A b = — S (0) = 
o A 



Bp, 

0 



(3.20a) 



R.4 (3.206) 

Equation (3.20a) is a set of Normal equations for a 1-D multichannel problem which we 
solve for the prediction matrices A s and the prediction error covariance E P[ . Equation 
(3.20b) represents Normal equations for a 1-D single channel problem that we solve for the 
filter parameters and the prediction error variance a A . The solution of these two sets 
of Normal equations then allows us to compute the filter coefficients as A s ® a A and the 
error covariance as 



E, — a \ Er 



(3.21) 



In the case where the covariance matrix in (3.12) is completely separable, along rows 

and columns and between channels, we have R = R c ® R B ® R^ . Then we are led to the 
two 1-D single channel subproblems 



R>i i 



R B a J3 — °~B L 



B 



and the 2-D multichannel parameters are computed from 



A = I 



M 



® &A 



°l R C 



(3.22a) 

(3.226) 



(3.23a) 

(3.236) 
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The decomposition into lower-order linear prediction problems discussed here would 
seem to be restricted to the case where the multichannel 2-D random process is separable. 
In actuality, a decomposition of the multichannel 2-D linear prediction problem into lower- 
order problems exits in general where there is no separability. This is discussed in the next 
section. A difference arises however in that the lower-order subproblem for the general 
case involves an expanded multichannel problem. In addition, we are guaranteed in the 
separable case, that if the correlation matrix in (3.12) has doubly block Toeplitz structure, 
the correlation matrices of each of the subproblems will have Toeplitz or block Toeplitz 
structure. In the general case the Normal equations for the subproblems do not all involve 
correlation matrices with Teoplitz structure, even if the matrix in (3.12) has the required 
doubly block Toeplitz form. 

3.4 Relation Between Multidimensional and Multichannel 
Linear Prediction 



There is a close relation between single channel 2-D linear prediction problems and 
multichannel 1-D linear prediction problems (7j. More general results exist that show 
that multidimensional linear prediction problems of any dimension can be decomposed 
into a series of 1-D multichannel and single channel linear prediction problems [8]. Some 
specific results will be discussed here that show how 2-D multichannel problems relate to 
higher order 1-D multichannel problems. These results will be used in the next section 
to formulate a method for estimation of the 2-D multichannel model parameters without 
explicit prior estimation of the correlation function. 

Fig. 3.8 illustrates a multichannel 2-D linear prediction problem with first quadrant 
support (a) and the corresponding multichannel 1-D linear prediction problem (b). The 
Normal equations for the multichannel 2-D problem are given by Eqs. (3.12), (3.13) and 
(3.10). The structure of the equations is further detailed in Fig. 3.3. 

In Fig. 3.8(b) the data is regarded as an array of P 2 M channels of 1-D signals evolving 
along the n x direction. The term x n used in the index-ordered representation of the 2-D 
data corresponds to the signals in this array of channels. The 1-D P 2 M-channel linear 
prediction problem can then be written as 



Pi-i T 

( a(<) ) 2 C «-1 ( 3 - 24 ) 

,'= 1 

where x n is the P 2 -^f-diniensional data vector and x n is its estimate. The coefficients a 4 '* 
have the form 





a (i) 
w 0 0 


<*<•> 

U 01 


a' 0 l 

“o,P 3 -l 




oc ( 1 ) 1= 


a**) 

U 10 


. 


( n ) 


(3.25) 
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Figure 3. 




(a) 2-D dota 



1st prediction direct icfn 
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2nd prediction direction 
(b) Multichonnel data 



Multichannel 2-D linear prediction and related multichannel 1-D problem. 
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where each block a.['} is of size M x M. The a (n) are found by solving the equation 



(3.26) 



where E Pl is the P 2 M X P 2 M prediction error covariance and where the correlation matrix 
appearing on the left side of (3.26) is the same as that in (3.12). Now post mulitply both 
sides of Eq. (3.26) by the term Al 0 ' and compare the result to Eq. (3.12). In view of Eqs 
(3.10) and (3.13), Eq. (3.26) will be identical to Eq. (3.12) if we require 

E P ,A (0) =S (0) (3.27) 

and 

A (i) =a (i) A (0) ; i = 0, 1, . . . , — 1 (3.28) 

The foregoing equations show that the multichannel 2-D Normal equations can be 
solved by the following steps. First solve the P 2 M-dimensional 1-D multichannel problem 
(3.26). Since this problem is 1-D, the multichannel Levinson recursion can be employed 
to find the coefficient matrices and the prediction error covariance matrix E Pl . Next 
solve (3.27) for A*°E this is also a set of 1-D Normal equations although the matrix E Pl 
is not in general Block Teoplitz. Finally find the multichannel 2-D coefficients from (3.28). 

Since the multichannel 2-D coefficients are expressed as a product of terms in (3.28), 
we have the following interpretation of the multichannel 2-D linear prediction problem. 
The multichannel 2-D prediction error can be computed by first predicting the data along 
the direction shown in Fig. 3.8(b). This is a 1-D linear prediction problem and results in 
a P 2 M- dimensional prediction error vector corresponding to the array of P 2 M channels 
shown in the figure. This prediction error is itself filtered in the second direction to 
compute the M-dimensional prediction error vector for the 2-D problem. The Normal 
equations solved to obtain the filter coefficients for this second linear prediction problem 
are represented by (3.27). 

This view of multichannel 2-D linear prediction has one further interesting aspect. If 
the multichannel Levinson recursion is used to solve (3.26), both forward and_ backward 
1-D linear prediction parameters are generated as part of the recursion. By arguments 
similar to those presented above, the backward 1-D parameters can be shown to relate 
to the coefficients for the second quadrant 2-D filter. Thus an algorithm based on the 
multichannel Levinson recursion can be used to solve for both the first and second quadrant 
filter coefficients simultaneously. This will be shown to be of particular use for the spectrum 
estimation problem considered in Chapter IV. 
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3.5 A Direct Method for AR Parameter Estimation 



The connection between multidimensional and multichannel linear prediction suggests 
a method for estimating the AR model parameters. This method will be called a direct 
method since it does not require prior estimation of the correlation function. 

Our method capitalizes on the existance of a direct method for solving the 1-D mul- 
tichannel linear prediction problem. We will refer to the 1-D method as the multichannel 
Burg algorithm since it is based on ideas originally suggested by Burg [9]. Details of the 
method however were developed separately by Nuttal flOj and Strand [ 1 1] . The multi- 
channel Burg algorithm is described in Appendix D. 

In order to apply the multichannel Burg algorithm to 2-D data, the data is first 
partitioned into strips along the n 2 direction as shown in Fig. 3.9(a). The strips of width 
P 2 for an Pi x P 2 quadrant filter are catenated along the nj direction as shown in Fig. 
3.9(b). As in the previous section this data is considered to be a 1-D process (in the n x 
direction) with P 2 M channels. The multichannel Burg algorithm is then used to estimate 
forward and backward linear prediction parameters (including the error covariances). This 
procedure makes no explicit estimate for the correlation matrix. Discontinuities where the 
strips were catenated together are ignored since they represent only a small portion of 
the data. The error covariance matrices that are computed as part of the multichannel 
Burg procedure are used to form (3.27) and an analogous equation for the backward case. 
These equations are solved by conventional methods and the multichannel 2-D parameters 
are computed from (3.28) and its analog for the backward case. This completes the 2-D 
estimation procedure. 
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(a) 2-D data organization 
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Figure 3.9 Sectioning data for the direct method of AR parameter estimation 
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IV. Multichannel 2-D Spectrum Analysis 
4.1 Spectrum Estimation Models 



Parametric or model-based approaches to spectrum estimation are based on a 
model for the spectrum involving a finite and relatively small number of unknown 
parameters. These parameters can be estimated from the data and used in a formula 
that gives the spectrum of the model in terms of the parameters. 

Autoregressive modeling is one form of spectrum estimation that will be examined 
here. A model of the form of Fig. 3.7(b) is used to represent the data. Filter parameters 
and the white noise covariance are estimated by the techniques described in Chapter 
III. Since the input spectrum is assumed to be that of white noise, a formula based on 
Eq. (2.25) can be used to compute the spectrum. If 

H.(z,,z,) = A l.t, (4- 1 ) 

( 1 1 , t-z ) £ a 



and the output is assumed to be white noise with constant spectral matrix S v ^ , then 
the spectrum estimate is given by 



S (^! , u 2 ) — 1 (^! , u 2 ) H u T (a?! , lo 2 ) 

= H z 1 (z x j z 2 ) E w Hj (a x 1 , z 2 1 ) 



zi = e 






(4.2) 






Note that this results in an estimation of the entire spectrum matrix i.e., all of the 
auto spectra and cross spectra at once (see Eq. (2.2l)and Appendix A). If desired, 
normalized quantities such as the magnitude squared coherence can be computed from 
the elements of the spectral matrix (see Eq. (2.22) and Appendix A). 



Various spectrum estimates can be obtained by assuming different support for the 
AR model. These are discussed below. 



4.1.1 Non-svmmetric Half Plane Models 

Non-symmetric half plane models are a. suitable choice for the spectrum estima- 
tion since arbitrary 2-D spectra can be factored into sections with infinite extent non- 
symmetric half plane support [5]. Although any practical algorithm cannot use infinite 
or even very large support regions, nonsymmetric half plane models of moderate size 
have proven to give reasonable results for spectrum estimation. An interesting fact is 
that it is unnecessary to consider the forward and backward predictors of Fig. 3.5(a) 
and (b) spearately. While these filters have different parameters, the white noise covari- 
ance is related such that formula (4.2) gives theoretically identical spectrum estimates. 
In subsequent sections of this chapter we examine results of NSHP spectral models of 
various orders. The definition of order used for NSIIP models is depicted in Fig. 4.1. 
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1 Definition of order for NSHP models . 
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4.1.2 Quadrant Models 



AR models with first or second quadrant support are of interest primarily because 
of their convenience of computation and estimation of the model parameters. The third 
and fourth quadrant filtes, are also distinct but when they are used in (4.2) with their 
corresponding white noise covariance they give estimates identical to those of the first 
and second quadrant filters. A method for estimating the model parameters directly 
from the data without prior explicit estimation of the covariance function was described 
in Section 3.5. 

It will be seen that first and second quadrant AR models when used individu- 
ally give poor estimates of the spectrum of many random processes. Features of the 
spectrum such as peaks tend to be displaced in frequency and elongated in one direc- 
tion. However a certain combination of the models has proven to give good results 
for estimation of spectra in 2-D single channel problems ( Jackson & Chien [12] ]. We 
propose here a generalization of this method to 2-D multichannel spectrum estimation. 
A combined spectrum estimate is computed from 

S(wi,w 2 ) = 2 + S, - / (w : ,w 2 )) 1 (4.3) 

where S t and Su are the spectra corresponding to AR models with first and second 
quadrant support. From Eqs. (4.3) and (4.2) one can therefore write 

S [oji ,u>2) — 2 (^Hj (uq , u>2 )D lv ) Hj (uq , w 2 ) + H n (wi,w 2 )Ti w 1 jj Hj j [u>i ,u>2)^j (4-4) 

In the special case that the noise is normalized so that Ylw, = Ylw ,, = ^ Eq. (4.4) 
has the simpler form 

5(wi,w 3 ) = 2 {oj 1 ,u 2 )H i (oq ,w 2 ) + (uq , w 2 )H U (w x , w 2 )) (4.5) 

This form is analagous to the form proposed by Jackson and Chien for the single channel 
case. 

4.2 Resolution and Phase Estimation Experiments 
(Sinusoids in Noise Background) 

Simulated multichannel 2-D random signals consisting of sinusoids with various 
relative phases and frequencies in noise backgrounds were generated. This section 
presents estimated spectra of these random fields. 

4.2.1 Comparison of NSHP and Quadrant Modeling 

Here we discuss the estimation of a sinusoid in noise background and compare 
the results of NSHP modeling to quadrant modeling. Three numerical examples are 
presented which illustrate the results. 
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Example 1 



This example is concerned with the analysis of spectra for a two-channel 
2-D single sinusoid signal with different phase in additive noise. The signals 
generated in channels 1 and 2 were 

x l (n x , n 2 ) = cos(n 1 w 1 + n 2 w 2 ) +Wi (% , n 2 ) 

x 2 (n 1 , n 2 ) = cos(n 1 u> 3 + n 2 w 4 +(/>) + w 2 (n 1 , n 2 ) 

where w 1 (n 1 ,n 2 ) and w 2 (n, , n 2 ) are zero-mean independent white noise sig- 
nals. Spectrum estimation results are given for a data set size of 64 x 64, 
and a third order NSHP filter. Two cases are considered in this example. In 
the first case we suppose that, the two channels have the same frequency, but 
different phase 



'^i _ y, = -,u 3 = uq = -<p = 1 radian j , 

while in the second case, we assume two channels with different frequency and 
different phase: 



= u> 2 = — ,u 3 = u >4 = — and <f> = 1 rad^ , 



Fig. 4.2 shows the results for the components of the 2x2 spectrum matrix in the 
first case. Only the cross term S 12 (u> 1 ,ui 2 ) is shown (magnitude and phase) since the 
term S 21 (u 1 ,u 2 ) is theoretically and numerically identical. The results show a distinct 
peak in each of the three components S l± , S 12 , S 22 corresponding to the location of 
the sinusoid. The center of the peak is accurately located near (tt /2, tt/ 2). Although 
the phase of the cross spectrum shows various artifacts around the edge of the region 
(where the magnitude is small and the phase is that of the noise) the phase at the 
location of the sinusoid is accurately estimated. 

Fig. 4.3 shows 5'(w 1 ,w 2 ) for the second case (different frequency). The power 
spectrum estimates of the first and second channels *S' 11 (w 1 ,w 2 ) as S 22 (u 1 ,u 2 ) have a 
single-peak at the location of the sinusoid. Although the cross spectra should theoreti- 
cally show no presense of sinusoids some small amount of energy is detectable at those 
frequencies in the cross spectrum. Similar effects have been observed in 1-D multichan- 
nel spectrum estimation and have been attributed to non exact pole-zero cancellations 
in the estimate for the cross spectrum [13]. 
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Figure 4.3 Estimate of spectra for sinusoids of different frequency 
(NSHP model). 



42 - 



Example 2 : 

Two-channel 2-D signals with two sinusoids in additive noise are considered 
in this example. The signals in the channels are defined by 

x i (nj , ra 2 ) = cos^u;, — n 2 w 2 ) + cos(n 1 u> 3 4- n. 2 w 4 ) + w l [n l , n 2 ) 
x 2 (n l ,n 2 ) = cos(n 1 w 1 + rc 2 u > 2 + <f>i ) + cos(niU ; 3 + n 2 u> 4 + 4> 2 ) 

+ w 2 (n.i , n 2 ) 



The results are presented here for the following data: u», (n x , n 2 ) and i v 2 [n l ,n 2 ) 
are independent zero-mean white noise signals, the frequencies and phases are 
u>i = u> 2 = = cp 4 = j, and <f>i = cf> 2 = 1 radian and the data set size is 

64 x 64. 

The power spectrum estimation results for a second order NSHP model are given 
in Fig. 4.4. The results are close to the theoretical results. S u (wi,o/ 2 ) and S 22 (uq , u 2 ) 
shows that the two sinusoids are easily resolved. The estimated amplitudes are unequal, 
but this characteristic has been observed in even 1 - D AR spectrum estimates. The 
cross-spectrum S 12 (t<q ,u > 2 ) shows the sinusoids resolved and the phase estimate is close 
to the true phase of 1 radian. 

Example 3: 

Three sinuosids in each channel are considered in this example. The location 
of the sinusoids is shown in Fig. 4.5. Observe that the phase of the sinusoids 
in channel 2 differ from those in channel 1 by 1 radian. The sinusoids are 
imbedded in white noise as in the previous examples. Again a data set of size 
64 x 64 was used. 

Fig. 4.6 shows the results of spectrum estimation using a fourth order NSHP 
model. The results show good estimation of the position of the sinusoids in the au- 
tospectra and the cross spectra and there is almost no evidence of energy from sinusoids 
uj b and ui D appearing in the spectrum of the opposite channel or in the cross spectrum. 
The phase of the cross spectrum in the region where the sinusoids are located is nearly 
constant with a correct value of 1 radian. 

Spectrum estimates for Examples 1 and 2 were also computed using quadrant plane 
models. The estimates were based on 2 x 2 regions of support for the first and second 
quadrant filters. Fig. 4.7(a) shows the component 5 X1 of the spectral matrix. Results 
for S 22 are similar. The use of either the first of second quadrant filter alone results in a 
spreading of the peak in one direction. The combined estimate (Eq. (4.3)) gives a more 
accurate result similar to that of the NSHP model. Fig. 4.7(b) shows the magnitude 
and phase estimates for the cross spectrum S i2 . A similar spreading phenomenon is 
observed in the magnitude estimate but the estimation of phase is correct for both the 
individual and the combined spectrum estimates. 
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Fig. 4.8 shows the spectrum estimate for two sinusoids in white noise corresponding 
to Example 2. In this case the phases 4>i and <+> 2 were taken to be unequal and given 
values of 1.5 and 0.5 radians respectively. The results of this quadrant based model 
can be qualitatively compared to the results for the NSHP model of Fig. 4.4. The 
placement of the sinusoid along a diagonal shows some characteristics of the quadrant 
models. The first quadrant results show a very significant spreading of the peaks along 
a direction orthogonal to the line connecting their centers. This is observed in both 
the autospectral components and the cross spectra. The second quadrant estimates 
show good resolution with little spreading of the peak. However the combined estimate 
gives the best results with sharp peaks and with magnitudes more nearly equal than 
those observed with the NSHP model. The phase in Fig. 4.8(c) is slowly varying in the 
region of the sinusoids for all of the estimates with a correct values of approximately 
1.5 and 0.5 radians at the locations of the sinusoids. (Exact values produced by the 
combined estimate are 1.42 and 0.57 radians respectively.) 

These experiments indicate that the result of spectrum estimation using a single 
quadrant model is not generally reliable but the estimate resulting from combining the 
two models according to Eq. (4.3) is quite accurate. 

4.2.2 Model order. Data Set Size, and Signal-to-Noise 

Ratio Experiments 

A comprehensive set of experiments was performed to determine performance of 
the spectrum. estimation procedures as a function of model order, data set size, and 
signal-to-noise ratio. The results are briefly summarized here. More detailed descrip- 
tion of the results will be reported separately. 

The results of the model order and data set size experiments showed that to some 
extent the lack of resolution resulting from a small data set size could be compensated 
for by choosing a larger model order. There is a limit to this trade-off however since 
a larger model has more parameters and thus should require more data to estimate 
parameters that are statistically reliable. The experimental observation may be better 
restated in the following way. When the data set is large a smaller order model can 
produce results that are comparable to a large model. 

The case of two sinusoids in noise (Example 2) was repeated for smaller size data 
sets using a NSHP model. For a second order filter the results of spectrum estimation 
using a 32 x 32 point data set were essentially the same as those using the 64 x 64 
point data set reported above. For smaller data sets the results degraded considerably. 
However the resolution obtained using a second order NSHP filter on a 16 x 16 point 
data set was similar to that obtained with a third order filter on an '8 x 8 point data set. 
These results for the cross spectrum are shown in Fig. 4.9. Results for the autospectra 
are similar. Note in Fig. 4.9 that although the amplitude of the estimate varies in the 
three cases depicted, the phase remains essentially constant. 
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Figure 4.7 Estimation of single sinusoid using quadrant models. 
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Fipure 4.8 Estimation of two sinusoids using quadrant models. 
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Figure ^.8 Estimation of two sinusoids us .i tip; duadrant models. (eontM) 
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Figure ^.8 Estimation of two sinusoids using quadrant models. (corit'd) 
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The signal-to-noise ratio for sinusoids in noise was defined by 

SNR = 10 log 10 ( 4 - 6 ) 

where A is the amplitude of each sinusoidal component and o 2 is the white noise 
variance (chosen to be the same for each channel). Fig. 4.10 shows results of the cross 
spectrum estimation for a 64 x 64 data set with closely spaced sinusoids of frequencies 
uq = u> 2 — 7r/2 and uq = u> 2 = the model used for these experiments was second 
order. At a SNR of 12 dB the sinuosids are completely resolved with sharp peaks. At 
a SNR of 3.5 dB the peaks begin to merge and at 0 dB the peaks become a single 
ridge making it difficult to predict that there are two sinusoids. The phase obtained by 
this method however remains essentially constant for all of the signal-to-noise values. 
A contour plot of the phase is shown for the 0 dB case again illustrating its slowly 
varying character in the region around the two sinusoids. 

4.3 Estimation of More General Spectra 

This section contains two additional examples of multichannel 2-D spectrum anal- 
ysis involving simulated data. These cases were designed to test the accuracy of esti- 
mating a linear phase term in the cross spectrum and the ability to estimate parameters 
of a known 2-D random process from data. The latter is closely related to the 2-D lin- 
ear system identification problem, since a linear system can be identified by driving it 
with white noise and then estimating the cross spectrum between input and output. 

4.3.1 Estimation of Delay 

For this experiment two channels of data were defined by 
Zi (n!,n 2 ) = i ^(nj.na) 

x 2 (n 1 ,n 2 ) = pXi (n x - d 2 ) + w 2 [n l ,n 2 ) 

where W 1 and W 2 are two independent white noise processes. Spectrum estimates were 
generated and the slope of the phase of the cross spectrum was measured to estimate 
the delays d x and d 2 . Fig. 4.11 shows the cross spectrum estimate that was obtained 
using a first order NSHP model for white noise term with unit variance and parameter 
values of (3 = 0.5, d 1 — d 2 — 1. The magnitude of the cross spectrum is constant 
and the phase shows a linear dependence with slope corresponding to = d 2 = 1. 
Thus the experimental result agrees precisely with the theory. 

4.3.2 Estimation of Parameters for a Linear Model 

The goal of this experiment was to estimate the parameters of a known linear 
process from records of data generated by the process. The following multichannel 2-D 
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Figure 4.11 Linear phase estimation. 
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where: Wi(ni,n 2 ) (t = 1,2) are two independent white noise signals. This process 
actually has support only in the first quadrant but a first order NSHP model was used 
to test the estimation procedure. For this problem one would expect the two matrix 
coefficients A lt - i and A n to be small (near zero). 

Table 4.1 shows the result of estimating the parameters with a 64 x 64 point 
data set. The estimates for the non-zero parameters of the model are close to the 
given values, and all of the remaining parameters but one are at least one order of 
magnitude smaller. The estimated spectrum components for this process are shown in 
Fig. 4.12. Energy is spread over a wide range of frequencies with concentration near 
higher values of uq and lower values of w 2 . This is consistent with the signs of the 
terms in the defining equations for the process. 

Table 4.1 Estimated parameters for a 1st order NSHP 
model for a linear process. 
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4.4 Comparison of Direct and Indirect Methods 

The results described in the previous sections were based on algorithms that first 
estimate the 2-D matrix correlation function and then solve Normal equations to de- 
termine the AR model parameters. These methods will be called indirect since they 
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require estimation of the correlation function as a prerequisite to determining the model 
parameters. A direct method for estimating the model parameters, i.e., a method that 
estimates the parameters directly from the data was described in section 3.5. This 
method is very well suited to spectrum estimation using combined first and second 
quadrant models since it estimates both sets of filter parameters simultaneously. Al- 
though the direct method has not been exhaustively tested on all of the cases reported 
earlier, our initial results indicate that its estimates are at least comparable to and in 
some cases better than those of the indirect methods. 

Figure 4.13 shows the results of estimating a single sinusoid in white noise using 
combined 2x2 point quadrant filters. These results are the same as the combined 
results shown in Fig. 4.7 but are depicted in a slightly different format. Fig. 4.14 
shows the spectrum estimation results for the same data using the direct method of 
parameter estimation. The direct method yields a sharper peak with somewhat lower 
sidelobes. Estimates of the one radian phase shift between the channels are 1.05 rad. 
for the indirect method and 1.03 rad. for the direct method. 

Figure 4.15 shows the results for estimating two sinusoids in white noise using com- 
bined 3x3 point quadrant filters with parameters estimated by the indirect method. 
These results are the same as those shown earlier in Fig. 4.8. Fig. 4.16 shows the 
results of the direct method applied to the same data. In this case the results appear 
to be nearly equivalent. A few minor peaks appear in both cases. For the indirect 
method these peaks appear closer to the main peak while for the direct method they 
appear further out. Phase estimates at the peak locations are 1.47 and 0.55 radians 
for the indirect method and 1.48 and 0.50 radians for the direct method. The values 
produced by the direct method are slightly closer to the true values of 1.5 and 0.5 
radians. 
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Figure 4.14 Estimation of sinusoid in white noise; 
direct parameter estimation. 
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Figure 4.15 Estimation of two sinusoids in white noise; 

indirect parameter estimation. (cont'a) 
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V. Image Processing Applications 

This chapter deals briefly with application of the theory of multichannel 2-D signals 
to image processing. Three applications are discussed namely, segmentation of color 
images, image coding, and image spectrum analysis. 

5.1 Image Segmentation 

This section describes an algorithm for image segmentation that uses a multichan- 
nel 2-D model to represent the images involved. This type of model works best for 
images involving multiple textures, such as aerial photographs of the ground in rural 
areas. 

The application here will be to color or “false color” images (the latter are derived 
from images recorded on color infrared film). For these cases the multichannel 2-D 
signal has three components (M=3) corresponding to the red, green, and blue video 
signals (see Fig. 5.1). The same methods could be applied to image data from a 
multispectral scanning satellite that records data on four, seven, or even more IR 
channels. In this case the number of components M would be equal to the number of 
channels of scanner data employed. 

It is assumed that regions with similar texture are defined by boundaries of corre- 
sponding pixels in each of the three color image components and that within each region 
the texture is described by a multichannel 2-D AR model. A superimposed Markov 
model characterizes the occurance of regions by specifying a form for the probability 
that a given pixel belongs to the region given that neighboring pixels belong to the 
same or other regions. The combination of the two models has been called a doubly 
stochastic image model in analogy with doubly stochastic models occuring in the study 
of point processes. 

The doubly stochastic image model has been used in conjunction with monochrome 
images and single-channel AR models. Since many of the details of the algorithms are 
described in our earlier work [14] for monochrome images our description of the method 
here will be brief. 



5.1.1 Segmentation Algorithm 



The image to be segmented is assumed composed of several homogenous regions 
of texture. The red, green, and blue intensities are represented by the multichannel 
2-D signal* 



F(ni,n 2 ) = 



Fi{ni,n 2 ) 

F 2 (n L ,n 2 ) 
F 3 (n x , n 2 ) 



(5.1) 



* 



We use the notation F to represent an image to adhere to previous convention. 
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g. 5.1 Representation of a color image as a 3-channel 2-D random process. 
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and within a particular region this signal is generated by the AR model 

F' (n x , n 2 ) = Av.F'K - *1 ,n 2 ~ * 2 ) + W(n l5 n 2 ) (5.2a) 

* 1*2 

F (n l9 n 2 ) = F f (n l9 n 2 ) + p (5.26) 

where /x is a constant mean vector. The white noise W(n l5 n 2 ) is described by a 3 X 
3 spatially invariant covariance matrix 

= E [W(n x ,n 2 )W r (rci,rc 2 )] (5.3) 

which in general is not diagonal. A separate model of this type is formed for each of 
the image regions in which there is a different texture. 

The essence of the segmentation algorithm is as follows. Using the models (5.2) 
and a Gaussian white noise source, we can form a probability density function for the 
set of all pixels in the color image conditioned on the regions. Assume that there are 
Q regions R 1 , £ 23 • • - , Rq • Since the texture is independent from region to region, this 
density function can be written as 



Q 

p(F|^ 1 ,^ a ,...,^ 0 ) = JIp(F|^,.) (5.4) 

t= 1 

where p (F|£j) represents the joint density for the pixels contained in region £, . Now 
assume that there is a way to represent the probability of occurence of a particular set 
of regions in the image Pr [Z 1 , R 2 , . . . , Z Q \. Then Bayes’s rule states that 



Pr{R 1 ,R 2 ,...,R Q \F] 



p(F|£i£ 2 )---i^q) -Pr [#! , £ 2 , . . . Rq\ 

P(F) 



Thus a maximum a posteriori (MAP) estimate for the regions can be obtained if the 
Ri are chosen to maximize 



P { F|£ i5 £ 2 ,... R q ) Pr [Ri y R 2 ) - -.Rq] 



(5.5) 



Since the form of the probability density function for the white noise is known, 
the form of the density function for pixels within a given region can be obtained by 
transformation. The resulting image probability density function is most conveniently 
represented in terms of the errors of linear prediction which can be written as 

F(n.j , n 2 ) = F [n x , n 2 ) — ^ ^ t i2 F (n^ — , n 2 — z 2 ) (5-6) 

*i.*i 
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where the assumed constant mean level is subtracted out prior to applying (5.6). By 
using the linear model, the density function for pixels within a region can be repre- 
sented as 

p(F|£,)= Yl P» k , (E(ni,n 2 )) (5.7) 

where p Wk is the probability density function for the white noise source of type k { 
within region £, . 

Now, suppose that an image has many regions, but that each region contains only 
one or another of two texture types. Then choosing a set of regions for the image 
is equivalent to labeling the pixels with labels 0 and 1. If an appropriate statistical 
model for the labels of the pixels exists then the prior probability Pr , £ 2 , • • • £q] 
can be computed. This is the other ingredient needed to form the MAP estimate. An 
appropriate statistical model is one that represents a particular class of 2-D Markov 
processes [14, 15]. This model permits calculations of the desired prior probabilities. 

Combination of the linear filtering model with the Markov model results in an 
algorithm to obtain a MAP region estimate. Since the resulting equations are of high 
dimension and nonlinear, a (possibly suboptimal) solution is obtained by iterating the 
conditions 

[E<°»(n 1 ,n 2 )] T [E W0 ]- 1 E<°>(n 1 ,n 2 ) 

0 

+^|S Wo I - 2lnPr [0|S (ni , n ,)] > 

1 

E (1) (ni,n 2 )] T E (1) (n x ,n 2 ) 

+/n|E u , I |-2/nPr [l|S {rei ,„ 3) ] (5.8) 

which are derived from (5.4) through (5.7) and the Gaussian form of the density func- 
tion. The terms El°*(n!,n 2 ) and E (1 l(n 1? n 2 ) are the error terms at pixel (n l5 n 2 ) 
computed using the filters of type 0 and 1 respectively and Ea,. is the correspond- 
ing white noise covariance. The terms Pr [o| S ( „ x >ri 3 ) ] and Pr [l|5( raiin2 )] are Markov 
transition probabilities representing the probability that the pixel at location (n l5 n 2 ) 
has label 0 or 1 given that a set of other pixels S(„ i n3 ) in the neighborhood has a 
prespecified set of labels. The inequality provides a rule for assigning pixel labels indi- 
vidually based on labels obtained at a previous iteration; the two expressions in (5.8) 
are evaluated and the pixel is given a label corresponding to the sense of the inequality. 

5.1.2 Experimental Results 

Results of applying the segmentation to color images are given here. Fig. 5.2(a) 
shows a 128 X 128-pixel color image representing an aerial photograph of the ground 
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Figure 5.2 Color image of trees and field. 
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in a rural area.* The green textured area is a grove of trees and the predominantly 
yellow area is a field with tall grass. Figs. 5.2(b) and (c) and (d) show the red, green, 
and blue components of the image. (The blue component is very dark and may not 
reproduce well.) 

Figure 5.3 shows the results of segmentation of this image. The algorithm used 
3-channel 2x2 first quadrant models for each of the two terrain types. The Markov 
model used a neighborhood size of 5x5. Note that attempts to segment this image 
based on color would have failed because the fields contain large linear patterns of 
green. The linear filtering models however match the colored texture in the image and 
result in a segmentation, that except for a few isolated pixels, is very accurate. 

Figure 5.4(a) shows another color image of trees and a field. This image is some- 
what more difficult to segment because the areas of green in the field do not appear as 
a linear pattern but are mixed in more homogeneously. Figure 5.4(b) shows the results 
of segmentation using image models similar to those used in the other example. Again, 
with the exception of some small areas the segmentation is quite accurate. 

5.2 Image Coding 

The models developed in this report can be applied to the coding of color or other 
multichannel images for purposes of storage or transmission. In this application, a 
model is formed for the image, the parameters of the model are stored or transmitted, 
and the model is used to reproduce the image. Two general schemes appear to be 
practical. 

The first scheme is applicable when portions of an image each have a homogenous 
texture and one is concerned about reproducing only the general character of the texture 
and not the specific gray levels at each pixel. In this case a white noise-driven linear 
model can serve to characterize the image sufficiently so that only the parameters of 
the model need to be retained. This is analogous to the analysis-synthesis method for 
speech encoding [16] . In a practical application of the method, the image would first 
be segmented and the boundaries of the regions would be coded and retained. Within 
each region, a linear model could be formed for the image and the parameters (filter 
coefficients and noise covariance) would be estimated and retained. The image would 
be “decoded” by driving the filter for each region with white noise having the specified 
covariance and thus reproducing the image texture within that region. 

The second scheme is quite general and would apply to any image. It is analogous 
to differential pulse code modulation (DPCM) in speech. A linear predictive filter 
is derived and applied to the image. The coefficients for the filter and the actual 
error residuals are retained and coded. Since there is usually great redundancy in a 
multichannel image the error residuals will tend to have low dynamic range and can 

* Data courtesy of Rome Air Development Center, Griffiss AFB, N. Y. 
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Figure 5.3 



Segmentation of the tree-field color image 
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(b) segmentation 

Figure 5.4 Color image of another area of 

trees and field and segmentation. 
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therefore be coded with a fewer number of bits than the original image. The image is 
reproduced by inverse filtering using the decoded error residual as input to the filter. 

The method just described is most effective when the linear predictive filter is 
well matched to the image being coded. Since one cannot expect that large general 
images will be well-represented by a linear filtering model, it is advisable to divide the 
image into small rectangular blocks and develop separate linear predictive filters for 
each block. Since the image within each block is usually more homogenous, one can 
expect that lower error residuals will be produced, resulting in more efficient coding. 

DPCM methods have been developed and used in the coding of single images and 
image sequences. Their use for multichannel images using the linear predictive models 
discussed in this report would seem to be particularly effective. 

5.3 Image Spectrum Analysis 

Spectra for portions of the color image of Fig. 5.2 were computed using the 
methods of Chapter IV. Only the red and green components were analyzed; these 
were designated as channel 1 and and channel 2 respectively. 

Two portions of the image were selected: a 64 x 64-pixel section (one quarter of the 
image) in the lower left representing the field and a 64 x 64-pixel section in the upper 
right representing the trees. The signal index n x was taken as the horizontal direction 
and n 2 was taken as the vertical direction. The images showed distinct differences in 
their spectra which could be used as a basis for discrimination. 

Fig. 5.5 shows the 2-D spectral components of the fields. The red, green, and 
magnitude of the cross spectrum appear to be very similar. Power is distributed in the 
middle and higher frequency ranges in the uq direction and around zero frequency in 
the cj 2 direction. This power distribution can be attributed to the somewhat evenly 
spaced lines of green (perhaps some kind of plants) appearing in rows in the field. The 
phase of the cross spectrum seems to oscillate between approximately 10° and -40° in 
a more-or-less regular fashion. 

Fig. 5.6 shows the estimated spectral component for the trees. Again, the red, 
green, and magnitude of the cross spectrum are very similar. The power tends to 
concentrate in a lower region of the cj x scale than the power for the fields with a peak 
occuring around 7t/6 corresponding to the spatial placement of the trees. The phase 
of the cross spectrum again varies between +10° and -40° but with somewhat milder 
oscillations. 

The concentration of power around u; 2 = 0 in the tree spectra at first seemed 
strange. One would expect that the tree image would show similar intensity variations 
in the horizontal and vertical directions and so the spectra should be more-or-less 
similar in the cj x and u> 2 directions. On closer examination it was found that the 
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Spectra for red and green components of the fields image. 
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Fie 5.G Spectra for red and trreen comnoneiits of the trees irnaere. 



image intensity variations for this data are not the same in both horizontal and vertical 
directions. Fig. 5.7 shows typical slices of the image in the n x and n 2 directions (we 
have examined several). The intensity exhibits rapid changes in the nj (horizontal) 
direction but shows a very slow change in the n 2 (vertical) direction. This accounts for 
the energy in the spectra being concentrated around cu 3 = 0. 

We can only guess a reason for this lack of symmetry in the tree image data. A 
plausible explanation is that a slight blurring of the image occurred due to the motion 
of the aircraft during the time at which the photograph was taken. This would have 
resulted in a lack of higher frequency detail in the direction of motion. 
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Fig. 5.7 Spatial variation of intensity in the trees image. 



VI. Conclusions 



This report developed analysis methods for multichannel 2-D random signals. Al- 
though these signals could be considered as a special case of three-dimensional signals, 
the multichannel nature of their origin provides them with special properties that are 
not shared by multidimensional signals in general. A good example of this is the 
concept of phase of the cross spectrum. This idea seems natural when considering 
a multichannel 2-D signal but would have no counterpart in a 3-D analysis. While 
there has been a large body of research for two and three-dimensional signals, analysis 
methods specifically for multichannel 2-D signals have not heretofore been given much 
consideration. 

A substantial part of the work in this report concentrated on 2-D linear prediction 
and autoregressive modeling. Concepts from the analysis of both single channel 2-D 
signals and multichannel 1-D signals were generalized in this development. 

An important application of the signal modeling arises in spectrum analysis. The 
methods described in Chapter IV of this report showed howautoregressive models could 
be used to estimate the entire spectral matrix for multiple 2-D signals. Estimates are 
developed simultaneously for the autospectrum of each signal and the magnitude and 
phase of the cross spectra between the signals in each channel. Spectrum estimation for 
multichannel 2-D signals is the topic of a separate investigation [17] and many results 
in addition to those reported here have already been developed. 

The application of the analysis methods to image processing was discussed briefly. 
Methods for image segmentation, image coding, and spectrum analysis were discussed 
in this report. Many further image processing applications such as enhancement filter- 
ing, target detection [18], and others seem possible. 

While the report represents the results of a serious effort to develop and apply 
the theory of multiple 2-D random signals, the work is by no means complete. Several 
applications to image processing have already been mentioned and should be developed. 
In addition, the analysis of some types of radar and array data seem to be important 
applications of the spectrum analysis and modeling. Special sensors such as the laser 
radar [19] where the data collected represents intensity, range, and doppler in two- 
dimensions offer further opportunities for research. 

By concentrating on linear prediction and AR modeling, we have looked at only 
a single (although rich) aspect of the analysis of multichannel 2-D signals. Other 
more general methods of modeling and estimation including pole-zero and noncausal 
modeling were not fully developed. Representation of the signals in other ‘coordinate’ 
systems by applying transformations between channels may lead to data reduction or 
other advantages in various applications. We hope not only to continue our work in 
this area but also, by publication of this report, to stimulate interest in the further 
analysis of multichannel two-dimensional signals. 
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Appendix A 

Spectral Representation of Multichannel 2-D Random Processes 



A stationary multichannel 2-D random process can be characterized in the fre- 
quency domain by a 2-D spectral matrix. The spectral matrix is defined by Eq. (2.20) 
of Chapter II and has the specific form 



S x (c^i ,w 2 ) 



S x ! , to 2 ) 


s i a (wi,w 3 ) 


_ Sm i (^i ) ^2 ) 


Sm 2 (^ 1 , 072 ) 



$1M (^1 5 ^2 ) 

Sm M (tUl , w 2 ) 



(A.l) 



where S <7 - (u^ , w 2 ) is the 2-D Fourier transform of R i3 {k l ,k 2 ) the cross correlation be- 
tween channels i and j. The spectral matrix is a periodic function of uq and u 2 with 
period 27r in each dimension and so it needs only to be considered only on the interval 
— IT < (Jj 1 < 7T and — 7T < U) 2 < 7T. 



It follows from Eqs. (2.7a) and (2.11) that 

5,(wi,w a ) = Sf («„«,) 

and thus the spectral matrix is Hermitian symmetric. In addition, since the compo- 
nents Ra [k l , k 2 ) are positive definite functions, the diagonal terms S i{ (w! ,u> 2 ) which 
represent the 2-D power spectrum of each channel are non-negative. Cross spectral 
terms S i3 (o^ , oj 2 ) for i ^ j are complex in general. However, since R i} (Aq , k 2 ) is real 
the magnitude of S i} - is an even function and the phase is an odd function of and 
C^2 * 



It is sometimes convenient to write the matrix of the squared magnitudes of the 
terms of S T as 



Si i • • • o 



0 ... S 



MM J 



1 k\ 2 ... K, 



2 

1 M 



k m, ••• i 



[S 1X 



0 1 



M J 



where the middle quantity is called the squared coherency matrix. It’s elements are the 
magnitude squared coherency (MSC) between the i t/l and j t/l channels and are defined 
as 



«?i(w i,w 3 ) 



Su{ui,U2)S ]3 (u> i.,U> 2 ) 



(A.2) 



It follows from the properties discussed above that the matrix is real and symmetric, 
that is, (w! , u > 2 ) = ,w 2 ). The MSC measures the correlation between the the 

signals in channels i and j at frequencies ui 3 and u> 2 . It is insensitive to spatial shifts 
between the signals in the two channels. These effects are manifested as phase in the 
cross spectrum. 
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Appendix B 

Examples of Correlation Matrices 
This appendix gives examples of the correlation matrices 



R — E [x x T ] 



and 



R ' = E x'x' 



for index-ordered and component-ordered representations of a stationary multichannel 
2-D signal. The examples are given for Ni = 3,N 2 = 4, and M = 2. The terms r^ e 
appearing in the matrices are defined by 



r'l- = E [xi (n x , n 2 )xy (n x - k,n 2 - t)} 



That is, r’lj’ is the element in the i th row and the j th column of the matrix correlation 
function R(k,l) (see Eq. 2.6a). 
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Note : r.. = r.. 

'I Ji 

The matrix is ; Block Toeplitz but not block symmetric, has block Toeplitz 
(but not block symmetric) blocks. The sub-blocks are 
symmetric but not Toeplitz in general. 

Figure B.l Index-ordered correlation matrix = 3, N 2 = 4, M = 2* 
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Each major block is identical except for subscripts on r. The overall matrix 
is neither block Toeplitz nor block symmetric. Major blocks are block Toeplitz 
(but not block symmetric) with Toeplitz blocks. 

Figure B. 2 Component -ordered correlation matrix N ] = 3, N 2 = 4, 

M = 2. 
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Appendix C 

Direct Product of Matrices 



The direct product or Kronecker product of matrices is an important notational 
operation in dealing with multidimensional signals. This appendix summarizes the 
most important properties of the direct product without proof. For a more thorough 
discussion see Ref. [20]. 

Given two matrices 



On a 12 ••• cl 1 p 



A = 



and 



B = 



(CM) 



. a N i a N 2 



6n 6x: 



a N P 



f>.< 



(C.2) 



_b Ll b L2 ... b L Q _ 

the direct product B ® A is defined as the partitioned matrix 



E ? (g ) A . = 



Abi i 


Abi2 


Ab^Q 


Abi, x 


Ab L2 


. . Ab LQ 



(C.3) 



If matrices A, C and B, D are conformable then some properties are as follows. 



[B ® A) {D <g> C) = BD ® AC (C.4) 

If A and B are square matrices [N = P, L = Q) then 

(. B ® A)T = B T ®A t (<7.5) 

tr(B ® A) = tr{B)tr{A) (C'.O) 

\B ® A\ = \B\ n \A\ l (C. 7) 

Further if A and B are nonsingular then the direct product is nonsingular and 

{B ® A)' 1 = B' 1 ® A' 1 ' (C. 8) 



Some further useful properties follow from the preceeding ones. If A,$ and 0,^? 
are the eigenvalue and eigenvector matrices for B and A respectively, then A 00,$®^ 
are the eigenvalue and eigenvector matrices for B ® A. That is 

(B ® A) ($ ® <&) = ($ ® (A 0 0) (C.9) 
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To state this another way, if 4>i and are eigenvectors of B and A corresponding 
to eigenvalues A, and then 4>i 0 if),- is an eigenvector of B 0 A corresponding to 
eigenvalue A That is, 



(B ® A) ■ {fa 0 Tpj) - A i9j(6i 0 Vv) 



(C.10) 



In addition if A and B are symmetric with a modified 



Cholesky decomposition 



A = GDG t 



and 

B = G'D'G ,T 

where G and G' and unit lower triangular and D and D' are diagonal, then 



(B 0 .4) = (G' ® G)(D' 0 D)[G' ® G') T (C.lil 



Some further fundamental properties and cautions should be observed. It is clear 
from the definition that the direct product is associative, i.e., 

(C 0 B) 0 A = C 0 (B 0 A) (C. 12) 

In addition, the operation is distributive over addition: 

(B + C)®A = B®A + C®A (CU3) 

A 0 (5 4" C) = ^405 + 410^7 (C.14) 

However, the direct product is not commutative (in general B 0 A 40 B) and it is 
not distributive over ordinary matrix multiplication. 
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Appendix D 

Multichannel Levinson Recursion and Burg Algorithm 
Multichannel Levinson Recursion (l-D) 

The multichannel Levinson recursion (also called the Levinson-Wiggins-Robinson 
(LWR) algorithm) deals with a linear prediction problem of the form (3.24) and solves 
the Normal equations (3.26). The associated backward prediction problem is 



if 1-1 



x' = - V (a {i) 'Yx 

■X-n-Ni + l N x + l + i 



> = 1 



and leads to the Normal equations'* 



R ct'P 1 = R 



I 




p^i 


a (l)' 


= 


0 


_ CE 1 " ‘ - 1 > ' . 




0 



(3.24)' 



(3.26)' 



where a*’) has a form analogous to (3.25). The recursion is specified by the following 
set of equations. 



A, = [R r (l) R T (2) ... R T ({) ] d^_! (£>.la) 

a; =Af=[R(l) R(2) ... »(■)]«;„ (fl.lt) 

r< = (E;.,)"A, (£>•• 2a) 

r; =E 1 - , 1 Af (£>.26) 



a, -i " 




0 






0 




A-i. 

0 


r, 


(D.3a) 


0 






r; 


[D. 3b) 



e, = e,_, (i-r;r,) (dao) 

e; = b;_. (i-r,r;.) (£>.46) 

* Note that ^ denotes second level block reversal. 
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The terms F, and T( are the forward and backward reflection coefficient or partial cor- 
relation matrices. The term A, (A') can be interpreted as a cross correlation between 
the forward (backward) prediction error and the one unit delayed, backward (forward) 
prediction error. Eqs. (D.l) - (D.4) are applied recursively beginning with i — 1 and 
ending with i = P x — 1. 

Multichannel Burg Algorithm 

This multichannel form of the Burg algorithm was developed by Nuttal f 10] and 
Strand [ll]. The procedure centers around estimating the terms A, and A' of the 
multichannel Levinson recursion directly from the data. Then the rest of the equations 
(D.2) - (D.4) can be applied without modification. 

Define the forward and backward error terms for the 1-D multichannel linear pre- 
diction problem as 

§1*’ =x„ -x„ = ajx {D.5a) 

e'j 0 = x n _ ,•+ 1 -Z n -> + i = oc[ T x (D.56) 

From these definitions and Eqs. (D.3) it can be seen that the errors satisfy the recursion 

r r p '(*- i) 
dO 



P (0 - P (<-n 

Sn — in 



e '(0 - e '(>) _ p' e . 

—n + 1 a — n + 1 



‘ * —n 
.T 



(D.6a) 

(D.66) 



The recursion is sustained by defining 

ri= 1 

B = (e"^ 

n= 1 



[D.l a) 
[D.l b) 
(D.7c) 



where N P is the number of points in the data record minus the length of the filter. 
The matrix A { is then obtained as the solution of the bilinear equation 



where 



+ a,c 2 = c 3 


(£>.8) 


= (e;) _1 b 


(D. 9a) 


= (E 0“ l D 


(£>.96) 


= — 2G 


(£>.9c) 
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There are several approaches to solving (D.8). One approach is to transform A, into 
an equivalent vector representation by concatenating rows of the matrix into the MN X - 
dimensional vector 6. The matrix C 3 is likewise transformed to a vector c 3 . Equation 
(D.8) can then be written as 

Q6 = c 3 (D.10) 

where 

Q = C x <g> I w + I M <g> C 2 (D.ll) 

where 1^ is the M X M identity matrix. Eq. (D.10) can then be solved for 6. 



- 90 - 



References 



[1] R.M. Haralick and G.L. Kelly, “Pattern Recognition with Measurement Space 
Spatial Clustering for Multiple Images,” Proceedings of the IEEE, Vol. 57, No. 4, 
April 1969. 

[2] B.R. Hunt and O. Kubler, “Karhunen-Loeve Multispectral Image Restoration, 

Part I,” IEEE Trans. Acoustics, Speech, and Signal Processing , Vol. ASSP-32, 

No. 3, June 1984. 

[3] D.E. Dudgeon and R.M. Mersereau, Multidimensional Digital Signal Processing 
(Prentice-Hall, Englewood Cliffs, NJ, 1984). 

[4] A. Papoulis, Probablity, Random Variables and Stochastic Processes (McGraw-Hill, 

New York, 1965). 

[5] M.P. Ekstrom and J.W. Woods, “Two-Dimensional Spectral Factorization with 
Applications in Recursive Digital Filtering,” IEEE Trans. Acoustics, Speech, and 
Signal Processing, Vol. ASSP-24, No. 2, pp. 115-128, April 1976. 

16] T.L. Marzetta, “A Linear Prediction Approach to Two- Dimensional Spectral Fac- 
torization and Spectral Estimation,” Ph.D. Thesis, Dept, of Electrical Engineering 
and Computer Science, M.I.T. (Feb. 1978). 

[7] C.W. Therrien, “Relations Between 2-D and Multichannel Linear Prediction,” 
IEEE Trans. Acoustics, Speech, and Signal Processing, Vol. ASSP-29, No. 3, 
June 1981. 

[8] C.W. Therrien, “Linear prediction for 2-D Vector Random Fields,” Proc. 19th 
Asilomar Conference on Circuits, Systems, and Computers, Pacific Grove, CA., 
November 1985. 

[9] J.P. Burg, “Maximum Entropy Spectral Analysis,” Ph.D. Thesis, Dept, of Geo- 
physics, Stanford University (May 1975). • 

[10] A.H. Nuttall, “Multivariate Linear Predictive Spectral Analysis Employing Weighted 
Forward and Backward Averaging: A Generalization of Burg’s Algorithm,” Naval 
Underwater System Center, New London, CT, NUSC Tech. Rep. 5501, Oct 1976. 

[11] O.N. Strand, “Multichannel Complex Maximum Entropy (Autoregressive) Spec- 
tral Analysis,” IEEE Trans, on Automatic Control, Vol. AC-22, No. 4, pp. 
634-640, Aug. 1977. 

[12] L.B. Jackson and H.C. Chien, “Frequency and Bearing Estimation by Two-Dimensional 
Linear Prediction,” Proceedings 1st Int. Conf. on Accoustics, Speech and Signal 
Processing, Washington, DC, April 2-4, 1979, pp. 665-668. 

[13] S.L. Marple and A.H. Nuttal, “Experimental Comparison of Three Multichannel 
Linear Prediction Spectral Estimators,” IEE Proc., Vol. 130, No. 3 (April 1983). 

[14] C.W. Therrien, “An Estimation-Theoretic Approach to Terrain Image Segmenta- 
tion,” Computer Vision, Graphics and Image Processing, Vol. 22, No. 3, June 
1983. 

[15] J. Besag, “Spatial Interaction and Statistical Analysis of Lattice Systems,” Jo. 
Royal Statistical Society Ser. B, Vol. 36, No. 2, pp. 192-236, 1974. 

[16] M.L. Honig and D.G. Messerschmitt, Adaptive Filters: Structures, Algorithms, 
and Applications, (Kluwer, Boston, 1984). 



- 91 - 



[17] H. El Shaer, Dept, of Electrical and Computer Engineering, Naval Postgraduate 
School, Private Communication. 

[18] C.W. Therrien, T.F. Quatieri, and D.E. Dudgeon, “Statistical Signal Models and 
Algorithms for Image Analysis,” Proc. IEEE , Vol. 74, No. 4, pp 532-551, April, 
1986. 

[19] R.C. Harney and R.J. Hull, “Compact Infrared Radar Technology,” Proc. Society 
of Photo-Optical Instrumentation Engineers, Vol. 227 - C0 2 Laser Devices and 
Applications, pp. 162- 168, Washington, D.C., April 1980. 

[20] C.W. Therrien and K. Fukunaga “Properties of Separable Covariance Matrices 
and Their Associated Gaussian Random Processes, ” IEEE Trans. Pattern Anal, 
and Machine Intelligence, Vol. PAMI-6, No. 5, September 1984. 



- 92 - 



Principal Symbols 



x(n l5 n 2 ) 


Multichannel 2-D signal 


Zm (n l5 n 2 ) 


Signal component 


x,x n 


Signal vectors - index ordering 


5 


Signal vectors - component ordering 


x,x m 


Signal matrices 


A T . ,BT. 


Coefficients of difference equation, filter coefficients 


H{1 !,/ 9 ) 


Impulse response 


(*1 , ^2 ) 


System function (z transform of impulse) response 


(Wl.Wi) 


Frequency response (Fourier transform of 
impulse response) 


H, HK) 


Matrix representation of impulse response 


M 


number of channels 


m 


Mean for signal x 


R x (fcj , k,2 ) , i? T y (£]. , k 2 ) 


Correlation functions 


K x {k 1 ,k 2 ),K xy (k 1 ,k 2 ) 


Covariance functions 


S x (u 1 ,u 2 ),S xv (u 1 ,w 2 ) 


Power density spectra 


K x (Wi.Uj) 


Coherence function 


m 


Mean vector , index ordering 


m' 


Mean vector, component ordering 


R, R(n) 


Correlation matrix index ordering 


R',R ninj ,R n k ' n ' 


Correlation matrix, 
component ordering 


2, 


Multichannel 2-D prediction error covariance matrix 
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